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Abstract. - Kinetic theory of dissipative particle dynamics is developed in terms of a Boltz- 
mann pair collision theory. The kinetic transport coefficients are computed from explicit collision 
integrals and compared favourably with detailed simulations. Previous theory is found to 
correspond to a weak scattering limit, or Vlasov theory, and previously reported discrepancies 
with simulations are thereby resolved. In the large dissipation limit, we find qualitatively new 
scaling properties for the transport coefficients. 



Dissipative particle dynamics (DPD) is emerging as an attractive simulation method in 
such diverse areas as colloidal suspension rheology, explicit multiphase flow problems, polymer 
solution dynamics, and phase behaviour of block copolymer melts . The basic DPD fluid 
consists of a large number of identical particles interacting by pairwise soft repulsive, dissipative 
and random forces H . The soft repulsions typically correspond to an interparticle potential, 
U = A(l — r/r c ) 2 /2, where r is the distance between a pair of particles, r c sets the range 
and A sets the amplitude. The dissipative forces are similarly specified by f = —m^fWY)(r)v^ 
where m is the particle mass, 7 (with units of inverse time) sets the overall dissipation rate, 
the weight function is usually wu{r) = (1 — f/r c ) 2 , and v" is the relative velocity of the 
particles projected onto the line joining their centres. Given the dissipative forces, the random 
forces are completely determined by a fluctuation-dissipation theorem apart from an overall 
amplitude which determines the temperature k^T; their detailed form need not be further 
specified || || . All forces vanish for r > r c , and are arranged to conserve momentum locally. 
Thus hydrodynamics is recovered at longer length and time scales just as in a molecular fluid. 

For the applications cited above, the basic interactions are typically augmented by extra 
bond constraints, different types of particles, and so on. In nearly all cases though the 
properties of the basic DPD fluid are required to calibrate the results against other methods 
and to provide an essential bridge back to the real world. Thus there has been considerable 
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effort in establishing a rigorous theoretical basis for the method |3|] and in developing kinetic 
theories for the transport properties [js], |(|. These studies also open up relatively unexplored 
terrain for kinetic theory. For instance the 'ideal dissipative fluid' (basic DPD fluid with 
A = 0) is probably the simplest example of a completely structureless fluid with non-trivial 
transport properties. Dimensional analysis indicates that the transport properties of this fluid 
are completely determined by a dimensionless dissipation rate, 'yr c (k-QT/m) 1 ^ 2 , and number 
density, pr 3 , where p is the number of particles per unit volume. In what follows, we shall 
use units in which m — r c = k^,T = 1, so that the rms velocity of particles is v rms = \/i for 
instance. In these units, the general DPD fluid is completely specified by 7, p and A. 

Our starting point is the kinetic theory developed by Marsh, Backx and Ernst (MBE) 
||. They identify the principal time scales in the ideal dissipative fluid, and provide explicit 
expressions for the transport coefficients. For instance, the (normalised) MBE velocity auto- 
correlation function (VACF) is predicted to decay exponentially, 4>(t) = e~*/ r , with a decay 



For later use we define Ambe = ™ rms = 45\/3/27r7p to be a representative mean free path in 
the MBE theory. The self diffusion coefficient in the MBE theory, given our choice of units 
above, is simply D = r. The viscosity in the MBE theory has two contributions: a kinetic 
contribution, t]k — pD/2, and a dissipative contribution arising directly from the dissipative 
forces, 7]b = (7/9 2 /2) Jd 3 rr 2 wuir) = 27T7p 2 /1575. In simple terms these correspond to the 
internal friction induced by particles diffusing across streamlines, and the drag force that acts 
between particles on different streamlines. 

The MBE theory captures the basic kinetic phenomena in the DPD fluid, and is highly 
successful in explaining such apparent paradoxes as a decrease in viscosity for certain parameter 
ranges, as the dissipation rate 7 is increased. When compared in detail with simulations 
though, certain discrepancies were already noted by the original investigators (6), and were 
explored more recently by Pagonabarraga et al JtJ . Neither set of authors identify the origin 
of the discrepancies, although correlation effects are an obvious candidate. 

To investigate these discrepancies in detail, we carried out a systematic study of the ideal 
dissipative fluid by simulation ||. Typical results for the VACF are shown in Fig. 1. We find 
that 4>{t) is remarkably insensitive to p but depends significantly on 7, once a basic scaling 
prefactor pry is removed from the time dependence. For 7^5 the VACF decays as a single 
exponential, to an excellent approximation, over two decades of magnitude (after which the 
signal becomes swamped by noise). The decay rate is plotted against 7 in the inset to Fig. 1(a); 
it shows a systematic deviation from the MBE prediction. For 7 > 5 significant deviations 
from single exponential decay start to appear, indicating the early onset of correlation effects. 
Of course one would expect correlation effects to be present to some extent for all 7, since 
they lead for instance to the celebrated long time tail <p(t) ~ i -3 / 2 ||. 

We also studied the two contributions to the viscosity. Typical results are shown in Fig. 2. 
Although ryD appears to be well captured by the theory, t/k is again systematically wrong, by 
as much as a factor of three (Fig. 2(b)). It too is remarkably insensitive to density (Fig. 2(a)). 

When the VACF is well approximated by a single exponential, significant discrepancies 
with theory cannot be ascribed to the correlation effects hypothesised earlier. Furthermore, 
the lack of dependence of the kinetic properties on density is startling. For a long time these 
results puzzled us, until we recalled the basic properties of the classic Boltzmann pair collision 
theory ]To| . Let us first introduce an estimate of the number of collisions a particle undergoes 
before it loses its memory of its initial velocity. To be specific, define n co n = pX where A is the 
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mean free path (note that the collision cross section ~ r 2 = 1). In the MBE theory, this gives 
"■coil = /oAmbe — 45V3/27T7; we will use this estimate to analyse the data. The importance of 
a parameter proportional to /cAmbe was recently noted by Evans in a different approach . 

Our central argument is that n co n and pair collisions in general are the key to understanding 
the results. Firstly note that in a pair collision theory, the transport coefficients pD and 
?7k have no dependence on density. Secondly, when the VACF decay rate (Fig. 1) and t/k 
normalised by the MBE prediction (Fig. 2(b)) are examined as functions of n co u ~ I/7, a 
data collapse is found. Thirdly, both the VACF decay rate and ijk systematically asymptote 
towards the MBE theory in the limit n co n — > 00. These observations imply that the kinetic 
properties of the fluid are principally determined by pair collisions (n co ii in effect), and that 
MBE theory is actually only weak deflection theory, valid in the limit where many collisions 
are required to decorrelate a particle's velocity. 

If pair collisions are the determining factor, the transport properties ought to correspond 
to certain collision integrals in a Boltzmann theory. It is straightforward to show that the 
standard expressions hold for the DPD fluid with the minor modification that, for a given 
impact parameter and pre-collisional set of velocities, the random force leads to a distribution 
of outgoing velocities rather than a single, unique set as is the case for normal fluids. The 
collision integrals take the generic form M 



where b is the impact parameter in units of r c , v 12 is the pre-collisional relative velocity, v' 12 is 
the post-collisional relative velocity and 4>o is the Maxwellian distribution function for relative 
velocities. The functions are Fd(v,v') = v x (v' x — v x ) for the self diffusion coefficient, and 
F v (v, v') = v x Vy(v' x v' y — v x v y ) for the kinetic contribution to viscosity. The angle brackets 
denote an average over the distribution of out-going velocities for a given b and V12 . 

To evaluate these collision integrals, we used a numerical, Monte Carlo approach. A value of 
b 2 was chosen from a uniform distribution in the range [0, 1] and a value of V12 was taken from 
the Maxwellian distribution. We then followed the particle trajectory using the algorithm 
of Pagonabarraga et al [j7|. We noted the outgoing relative velocity and thus formed the 
integrand for this particular trajectory. By averaging over many trajectories we obtained the 
overall integral. Typically we used a time step of 0.005 to compute the trajectory and I0 8 
trajectories to get an precision of better than 1%. 

Once the collision integrals were determined, we calculated the diffusion constant and the 
viscosity using the standard results (in our chosen set of units) pD = — ^ and ?/k = — 
These results arc at a first Sonine polynomial level of approximation. For the hard sphere fluid, 
for example, this leads to an error of < 2%, compared to the true Boltzmann value ]I(J. To 
check whether the first Sonine polynomial approximation was as good for the DPD fluid, we 
also calculated the second Sonine polynomial correction. Just as in the hard sphere fluid, the 
correction was of the order of 1% so we believe the results quoted are excellent approximations 
to the true Boltzmann values. 

The predictions of the Boltzmann pair collision theory turn out to be in excellent agreement 
with the simulation results, thus confirming our premise that collisions are the key. Results 
for the decay rate of the VACF (defined by t = D) are shown as open circles in the inset to 
Fig. 1(b). Similar results for t]k in shown in Fig. 2(b). Moreover, we have been able to prove 
that the MBE theory is obtained in the Vlasov limit [flo| , where each collision only weakly 
perturbs the velocities of the colliding particles, thus verifying that the MBE theory is indeed 
a weak deflection theory. 

Now consider the significance of these ideas for the strong deflection limit. Since n co n ~ I/7 
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in the MBE theory, it appears that for 7 sufficiently large, n co \\ < 1. This is nonsense though, 
because at least one collision is needed to decorrelate a particle's velocity no matter how large 
7 is. Thus for 7 sufficiently large there must be a qualitative change in the behaviour. The 
crossover occurs at n co \\ ~ f or 7 » 7 c = 45\/3/27r w 12.4. Turning the argument around, if 
n C ott ^ 1 is enforced for all 7 there must be a change in the scaling of the mean free path from 
the MBE result, Ambe ~ 7c/(tp) f° r 7 ^ 7c to a new behaviour, A ~ l/p for 7 > 7 C . 

This argument implies that t, 13 and tjk should reach a plateau for 7 ^ 7 C - This is 
in marked contrast to the predictions of both the MBE theory and Evans' theory 
That a plateau is indeed found is demonstrated for tjk in the inset to Fig. 2(b). We find 
a good fit is t]k = 0.37(1) + 3.60(5)/7 (the numbers in brackets are estimates of the error 
in the final digit). Similarly for the self diffusion coefficient (data not shown here) we find 
D = 0.21(1) + 2.52(2)/7- At present we do not have an explanation for the straight line 
fits, but we can use them pragmatically to define crossovers at 7 ~ 3.60(5)/0.37(l) = 9.7(3) 
for Tj K and 7 » 2.52(2)/0.21(l) = 12.0(6) for D. These are satisfyingly close to j c » 12.4. 
This conclusion has implications for Schmidt number, an important dimensionless material 
parameter defined by Sc = rj/pD. In the large 7 (fixed p) limit, the viscosity is dominated by 
tjd ~ p 2 7 but the above argument changes the scaling of D from (p7)~ 1 to p^ 1 . This suggests 
that Sc ultimately grows as p 2 7 and not (P7) 2 as thought previously B. 

A pair collision theory is strictly valid only when the mean free path is large compared 
to r c , when effects such as finite collision volumes and correlations are unimportant. The 
excellent comparison with simulation implies though that the transport coefficients are rather 
insensitive to these effects so long as A > 1. There must be a point though where these effects 
become important, and we now turn to a discussion of this, the large p limit. 

Firstly, note that from a kinetic theory point of view the large p limit is rather unusual. 
For hard spheres for instance there is a natural maximum density, but for the ideal dissipative 
fluid, p can be increased indefinitely In this limit each sphere is in continual interaction with 
many others simultaneously and a picture in terms of pair collisions is surely invalid. Moreover, 
one can derive some of the MBE results by assuming a mean field interaction of this type . 
Therefore at very high densities, one would expect the MBE results (for r and D at least) to 
be recovered. A slightly more sophisticated approach would couple the motion of an individual 
particle to the hydrodynamic modes of the fluid, and recover for instance the celebrated long 
time tail in the VACF. An elegant mode-coupling theory along these lines has recently been 
developed by Espanol and coworkers |Q . 

Fig. 2(a) shows how 77K starts to drift down from the Boltzmann pair collision theory toward 
the MBE result for densities p > 10. It is at first sight remarkable that a pair collision theory 
works at such high densities. We can partially understand this observation though when we 
consider the fluctuations in t _1 in Eq. ([!]), interpreting the integral as a drag-per-particle in a 
random stationary background of other particles. If we suppose that the number of particles 
in a volume element d 3 r is a Poisson distribution with mean p d 3 r and variance equal to the 
mean, then r -1 is a weighted sum of independent random variables whose mean is given by 
Eq. (|l|) and variance by 

^r^ = \Jpd^[l^{T)f = ^. (3) 

Thus the ratio of the standard deviation to the mean for r _1 is (p//°c) ^ 2 where p c = 45/7n. 
For p — 10, for instance, this ratio is about 45%, and to get the ratio < 10% requires p > 200. 
Thus unless /> » 1 a particle sees large fluctuations in its interaction with other particles 
(O(l) in essence), and it is perhaps not surprising that a pair collision theory is applicable 
even though the apparent density becomes rather large. Another way of stating this conclusion 
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is to say that rj? severely over-estimates the effective interaction volume. 

The above discussions have concerned the ideal dissipative fluid. Let us finish off by briefly 
discussing the properties of the non-ideal fluid. The main effect of the soft repulsion forces is 
to provide another collision mechanism which is most significant at low 7; in fact at 7 = 
one recovers a regular fluid of soft spheres. The effect is to remove the I/7 divergence in 
the kinetic coefficients as 7 — * 0. Combined with the plateau that is expected for 7 > 7 C , 
this makes the transport coefficients very insensitive to parameter variations. The collision 
integrals can also be evaluated for non-ideal interactions, and again we find excellent agreement 
with simulations. This is so even though the VACF in some instances develops a pronounced 
non-trivial structure. 

In summary therefore, we have argued that, for the typical parameters used in applications, 
the kinetic transport coefficients in DPD are best captured by a Boltzmann pair collision the- 
ory, even though the apparent density is rather high. The evidence is the excellent agreement 
with simulation results, and the fact that the number of collisions a particle experiences in a 
mean free path appears to be the principal parameter that determines the kinetic properties. 
The kinetic theory of Marsh et al || is recovered as a weak scattering or Vlasov theory, valid in 
the low 7 limit. Previously reported discrepancies with this theory || Q are thereby resolved. 
Wc further argue that a qualitative change in the scaling of kinetic transport coefficients has to 
occur at large 7 since at least one collision per mean free path is required to decorrelate a par- 
ticle's initial velocity. This conclusion, too, is supported by our simulation evidence. Finally, 
the present study is likely to have implications for the analysis of transport behaviour in the 
various generalisations of DPD |ll| . More details of our results and our parallel investigations 
of the non-ideal fluid will be the subject of a longer paper. 

We acknowledge useful discussions and correspondence with P. Espanol, M. H. J. Hagen, 
C. Lowe and I. PAGONABARRAGA. AJM would like to thank the British Council and the 
Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) for providing a travel grant 
which aided this work enormously. 
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Fig. 1. - Velocity autocorrelation function for ideal dissipative fluid: (a) at p — 3 for 7 = 1 (steepest 
curve), 2, 4, 8, 16, and (b) at 7 = 2 for p — 1, 2, 3, 4, 5, 6, showing data collapse (note that a basic 
scaling prefactor p7 has been extracted from the time dependence). The dashed line in (b) shows the 
universal decay predicted by Marsh et al |6J (MBE). Where single exponential decay occurs over two 
decades or more, the inset to (b) shows the aecay rate normalised by the MBE result, tmbe = 45/27rp7, 
(points with error bars), as a function of 7. The open circles in this plot are predictions from the 
Boltzmann pair collision theory discussed in the text. 




Fig. 2. - Viscosity of ideal dissipative fluid: (a) the kinetic and dissipative contributions, t/k and rjo, 
at 7 = 4.5 as a function of density, and (b) data collapse for all tjk data (p = 0.25-10, 7 = 2-20), 
normalised with the MBE result, t^k.mbe = 45/47r7, and plotted against pAmbe = 45V3/27T7. The 
dashed lines in (a) and (b) are the predictions of MBE theory. The long dashed line in (a) and the 
open circles in (b) are predictions from the Boltzmann pair collision theory in the text. The inset in 
(b) shows all the r/K data plotted against I/7 with a linear regression line. Note the non-zero intercept 
in this plot, which indicates a plateau in rjK for large 7. 



